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Abstract 



o 

J>^, When inclusions with extreme conductivity (insulator or perfect conductor) arc 

closely located, the gradient of the solution to the conductivity equation can be arbi- 
trarily large. And computation of the gradient is extremely challenging due to its na- 
ture of blow-up in a narrow region in between inclusions. In this paper we characterize 
explicitly the singular term of the solution when two circular inclusions with extreme 
conductivities are adjacent. Moreover, we show through numerical computations that 
the characterization of the singular term can be used efficiently for computation of 

■ the gradient in the presence adjacent inclusions. 

< 
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1 Introduction 

> 

Frequently in composites which consist of inclusions and background (the matrix), the 
(f-) ' inclusions are closely spaced, and it is quite important from a practical point of view 

to know whether the gradient of the potential can be arbitrarily large as the inclusions 
I/"") . get closer to each other. The gradient of the potential represents the stress in anti-plane 

elasticity and the electric field in the conductivity problem; see [6]. It is known that the 
gradient of the potential may blow up as the distance between the inclusions goes to zero 
and their material parameters (conductivities or stiffness) degenerate. 

Suppose that B\ and B2 are inclusions whose conductivity is k. We suppose that the 
conductivity of the background is 1 (k 7^ 1). Let e be the distance between B\ and B2 
and assume that e is small. The problem is to estimate |Vtt|, where u is the electrical 
potential, in terms of e when e tends to 0. 
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There have been important works on this problem. If k stays away from and oo, 
i.e., c\ < k < C2 for some positive constants c\ and C2, then it was proved by Bonnetier- 
Vogelius [TO] and Li-Vogelius [19] that |Vtt| remains bounded regardless of e. This result 
was extended to elliptic system by Li-Nirenberg [18]. It is worth emphasizing that the 
results in |19l [T8] are not only for two inclusions case but also for the case of arbitrary 
number of inclusions. 

On the other hand, if k is either (insulating) or oo (perfectly conducting), then Vu 
may blow up as e tends to 0. For two identical perfectly conducting circular inclusions it 
was shown in [7] (see also [T7] and [JO]) that the gradient in general becomes unbounded 
as e approaches zero and the blow-up rate is e -1 / 2 . In [31 S], a lower bound and an upper 
bound for the gradient has been obtained. These bounds are valid for all k including 
extreme values (k = and k = oo) and provide the precise dependence of Viz. on e, k 
and radii of disks. The blow-up of the gradient may or may not occur depending on the 
background potential. In [5], Ammari et al characterize those background potential which 
actually make the gradient blow up. In |22U23| . Yun showed that the blow-up rate is e -1 / 2 
for perfectly conducting and insulated inclusions of arbitrary shape in two dimensions. In 
three dimensions, Bao et al [8] proved that the blow-up rate for the perfectly conducting 
inclusions is | e log e| — 1 and extended the result to the case of multiple inclusions [9]. Lim- 
Yun [20j also found the same blow-up rate when inclusions are spheres. Their estimates 
explicitly reveal the dependence on the radii of the sphere. They also showed in [21] that 
if there is a small bump in between two inclusions in two dimensions, then the magnitude 
of the blow-up gets larger. 

The purpose of this paper is to characterize the singular term of the solution, i.e., 
to establish an asymptotic formula for the blow-up of the gradient when two circular 
inclusions get closer. We find the decomposition of the solution u to the conductivity 
equation as 

u = g + b (1.1) 

where Vg may blow up at the rate of e -1 / 2 while Vb stays bounded regardless of e, when 
B\ and B2 are disks and k is either 00 or 0. We actually obtain an explicit formula for 
the term g which gives a precise description of singular behavior of Vu. 

The characterization of the singular term of the solution finds a very good application 
in the computation of electrical fields. Computation of the electrical field in the presence of 
closely located inclusions with extreme (0 or 00) conductivities is known to be a extremely 
difficult problem because of the the blow-up phenomenon in a very narrow region between 
inclusions. Since the the gradient of the solution is arbitrarily large, we need very fine 
mesh to catch the large gradient in a narrow region. The results of this paper constitute 
a significant step toward overcoming this difficulty since the singular term g is explicit 
and computation of b requires only regular meshes. We present efficient methods to use 
the decomposition for the computation of the solution and some results of numerical 
computation using them. Numerical examples of this paper show that these methods 
work pretty well. 

This paper is organized as follows. In the next section we derive the decomposition 
(jl.ip for the perfect conductors in the free space. In section 3, we deal with the same 
problem in bounded domains. Section 4 is for the insulators. New numerical methods and 
results of computation are presented in the last section. 
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The result of this paper can be extended to perfect conductors of spherical shape in 
three dimensions. This result will be presented in a forthcoming paper. 



2 Free space problem-perfectly conducting case 

Let Bj = B(cj,rj), j = 1, 2, be the disk centered at c,- and of radius r,-, and 

[k onB l UB 2 , 

which represents the conductivity distribution: the conductivity of the inclusions is k 
(k 7^ 1) and that of the background is 1. The equation we consider is 

V • o-Vu = in M 2 , (2.2) 

which may be viewed as the conductivity equation or anti-plane elasticity equation. A 
condition at the infinity is prescribed by 

it(x) - i?(x) = 0(|x| -1 ) as|x|^oo, (2.3) 

where H is an entire harmonic function and represents the background potential. 

If k = oo, the equation (|2.2|) with the condition (|2.3|) is understood as the following 
problem: 

Ait = in R 2 \ By U B 2 , 

u\dBj = Aj (constant) j = 1,2, (2-4) 

it(x) — H(k) = Odxj" 1 ) as |x| — > oo. 

The constants Xj can be determined by the additional requirements 

du 

— ds = 0, j = l,2, (2.5) 

dBj dv + 



where v is the outward unit normal vector of IR 2 \ B\ U B 2 , i.e., directed inward of Bi. 
Here and throughout this paper, the notations |+ and |_ are for limits from outside and 
inside inclusions, respectively. 

Let Rj, j = 1,2, be the reflection with respect to dBj, i.e., 

RjW := r f X ~ C £ +Cj, j = l,2. (2.6) 

|X Cj\ 

It is easy to see that the combined reflections R\R 2 and R 2 R\ have unique fixed points, 
say pi and p2, respectively. Let 

fc(x) := ^ (log |x - pi| - log |x - p 3 |) . (2.7) 

The function h, which was first found in [20J, has a special property: it is the solution to 

' Ah = in E 2 \ 5iUB 2 , 

h\dB 4 = Cj (constant) j = 1,2, 

dh ( 2 -8) 
ds = (-iy, 3 = 1,2, 



dBj dv 



k /i(x) = 0(|x|- 1 ) as |x| 



oo. 
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The following formula was proved in [22^ 123] : let Ai and A2 be constants appearing in 
(l2T4"j) . then 

A 2 -Ai= f Hd u hds = H(p 2 )-H( Pl ). (2.9) 

JdB 1 UdB 2 

The following is the first main theorem of this paper. 

Theorem 2.1 Let n be the unit vector in the direction of p 2 — pi and let p be the middle 
point of the shortest line segment connecting dB\ and dB 2 . For a harmonic function H 
in M 2 , let u be the solution to (|2.4]) . Then, the solution u can be expressed as follows: 



ti(x) = ah(x) + 6(x) (2.10) 

where 

a= 47TTir2 

n + r 2 

and for any bounded set £1 containing B\ and B 2 there is a constant C independent of e 
such that 

||V6|U« (n \ (fllUB9)) < C. (2.12) 
The asymptotic formula Vu as e — > is then given by 

V«(x) - JS3.(„ ■ VH)(p) f yi^, - + 0(1). (2.13) 

n+r 2 \\x-pi\ |x-p 2 | / 

Let us make a few remarks on Theorem 12.11 before proving it. It is shown in [22^ [23] 
that the fixed points pi and p 2 are given by 



p 1 =(-v / 2 A /^^v / ^ + O(e),0) and p 2 = (^U^-Jl + O(e), 

if ci = (— ri — 1,0) and c 2 = (r 2 + |,0). If r\ and r 2 are bounded below by a positive 
constant ro, then there are positive constants C\ and C 2 depending only on ro such that 



C ij^i <|VM*)|<C 2V t!±^ (2.14) 

for all x on the the shortest line segment connecting dB\ and dB 2 , see [21]. Thus the 
blow-up rate of |Vn| is e -1 / 2 . It is also proved in the same paper that 



V Ve 

for all x. Thus, an optimal bound for V« in a bounded domain can be obtained from 
(|2.14p and f|2. 15|) in terms of r\, r 2 , e and (n • Vi7)(p). In view of the formula (|2.1ip of a 
(we call it the stress intensity factor), the blow-up does not occur if (n- V-ff)(p) = 0. This 
fact was already found in [5]. One can also show that if r\ and r 2 are 0(e), then there is 
a constant C independent of e such that 

|V/i(x)| < j (2.16) 

for all x. Thus (|2.13p means that in this case, no blow-up occurs: stays bounded. This 
finding is in agreement with that in [3]. 

We first prove the following proposition by modifying an argument of Bao et al [5]. 
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Proposition 2.2 Let 

6(x) = u(x) - I J 02 ;\ UD1 U(x), xgR^^UBj), 

For any bounded set Q\ containing B\ and B2 and O2 containing Q,\, there is a constant 
C independent of e swc/i i/iai 

ll^llc 1 (ni\(BiUB 2 )) ^ Cll-^lli 00 ^)- ( 2 -18) 



Proof. It can be easily seen that b is bounded. In fact, since b is harmonic in R 2 \B\ U i?2 
and 

b\dB 2 - MdBi = 0, 

we infer from the result in pQ (see also [5j) that 6 is bounded in fix- 

Since h(x) — > as |x| — > 00, by the maximum principle, h attains its maximum and 
minimum on <9£?2 and dB\, respectively. Thus, we have 

II /i IIl°°(R 2 \( j B 1 UB2)) - h \dB 2 ~ h\dB!, 

and hence 

||fe||L°°(n 2 \(-BiUB2)) ^ IMlL°°(n 2 \(BiUB 2 )) + \ u \dB 2 - u\q Bi \ ■ (2.19) 

Let Xj = u\qb j — 1)2, and assume that A2 > Ai without loss of generality. Since 
(u — H)(x) —7- as |x| —> 00, the maximum and minimum of u — H occur on dB\ U 'dB<i- 
So, we have 

(u - iJ)(x) < max (u - H) < A 2 + \\H\\ Lx , {n2) , 
0B1U0B2 

and 

^1 — \\H\\L°°(no) < m i n (u — H) < (u — H)(x.) 
for all x£i 2 \ (B\ U B2). Since min R 2\( BlU £ 2 )(-u — -ff) < 0, we have 

Ai < ||-ff||£<x.(n 2 )- 

Therefore, we have 

\\u - -ff||L°°(r22\(BiUB2)) - ^2 + ||-H1|.L°°(0. 2 ) 

< A 2 - Ai + 2||F|| it x>(n 2 ) < 4||iJ||i<x.(Q 2 ), 
where the last inequality comes from (|2.9p . Thus 

\\ u \\L° a (n 2 \(B 1 uB2)) < Il-ff||i°°(n2) + \\ u ~ ^IU°°(n2\(BiuB2)) - ^\\^\\L^(n 2 )- 
It then follows from (f2T9|) and (j2~9|) that 

IHlL°°(n 2 \(-BiUB 2 )) < 7 ll- ff IU° (n 2 )- (2-20) 

We now show that 

l|V6||L°°(f7i\(BiUB2)) ^ CII-^IU 00 ^)- (2-21) 



5 



For that purpose, we define the harmonic functions G+ and G_ as follows: 



AG± = 0, in n 2 \ Bi U B 2 , 

G± = ±||&|Uo°(n 2 \( Bl uB 2 )) on dn 2 , 

G± = b, on dB 1 U a^ 2 . 



Then, ±(G± - 6) > in fi 2 \ -Bi U B 2 and G± - 6 = on dB 1 U dB 2 . By Hopf's Lemma, 
we have 

d v G + <d v b<d v G- on dB 1 UdB 2 . (2.22) 

We introduce more harmonic functions G + i, G +2 , G-\ and G- 2 defined as follows: for 
i = l, 2, 

' AG ±i = 0, in n 2 \ B h 

< G±i = G± = ±||&||l°°(o 2 \( Bi ub 2 )) on d£l 2 , 
k G ±i =G± = b, on dBi. 

Since b\gB 1 = &|aB 2 = constant, we have 

G+i(x) > 6| SBlUSB2 . 

In particular, 

G+i(x) > 6(x) = G+(x) on U aB 2 - 
Since G +i |an 2 = G + |^ 2 , we have 

G +i - G + > in ft 2 \ Bx U B 2 - 

Since G+j — G+ = on <9-Bj, it follows from the Hopf's Lemma that 

d u G +i < d v G + on dB h i= 1, 2. (2.23) 
Similarly, one can show that 

d v G-i>d v G- on dBi, i = l,2. (2.24) 

Note that G±i/||6||^cx>(q 2 \ m 1 uB 2 )) * s a harmonic function in VL 2 \ B\ which is ±1 on 
<90 2 and also has a constant value between —1 and 1 on dB\, and that dist(l?i, d$l 2 ) > cq. 
Thus, G-ti/||6||^oo(f7 2 \( BlUB2 )) can be extended as a harmonic function into fi 2 \ B{c\,r) 
where ci is the center of B\ and f is strictly less than the radius of B\ independently of 
e. Then, we have from interior regularity estimates for elliptic equations and (|2.20p that 

\\dvG±i\\L<x>(dB 1 ) < C||^IU°°(n 2 \(SiUB2)) - 7G||-£f || L oo(q 2 ). 
It then follows from (I232]) . ([2331 and (I2T24"]) that 

\\dvb\\L°°(dBi) < C||-^IIl°o(q 2 ) 
for some constant G independent of e. Similarly one can show that 

ll<9i^llL°°(aB2) - C\\H\\ L oo(Q 2 y 
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Since b is constant on dB\ and 8B2, we get 

l|V&||z,°°( 9BlU a B2 ) < C||iI|| L oo(n 2 ). (2.25) 
The standard interior regularity estimate for harmonic functions shows that 

l|V6|| L oo( 9fil ) < C||H|| B cx.(Q 2 ). 

The maximum principle now yields (|2.2ip . and the proof is complete. □ 
We are now ready to prove Theorem 12.11 

Proof of Theorem 12. 11 After translation and rotation if necessary, we may assume that 
ci = (— r\ — §,0) and c 2 = (r 2 + §,0). Then p = (0,0) and n = (1,0). It is proved in 
1221 [23] that 



pi=(- V2, ^^V~e + O(e),0) and p 2 = (^2, ^^Je + 0(e), o) . 
V V ri + ro / V \/ n + ro / 



r\ + r 2 / V V r i + r 2 

Therefore, we get from (|2.9p 



u\ dB2 -u\ dBl =2V2d Xl H(0,0)J-^-^e + O(e) (2.26) 
as e — > 0. On the other hand, one can see that 



h\dBa ~ h\aB! = —f=- \ j 1 +r " 2 \/e + 0(e) 

V27T V T\Ti 



Therefore, we get from (|2.17p 



2V2d xl H{v)J^r 2 V~z + 0{e) 

U(X) = ; V /l(x) + 6(x) 

47rrir2 



n + r 2 



d Xl H{p)h{x) + 0(v^)/»(x) + 6(x). 



Note that the gradient of 0(y / e)/i(x) term is bounded because of (|2.15p and so is 6(x) by 
Proposition 12.21 Thus we obtain (|2.10p by setting 0(y / e)/i(x) + 6(x) to be the new b(x). 
This completes the proof. □ 



3 Boundary value problem- perfectly conducting case 

Let O be a bounded domain with C 2 -boundary containing two circular perfectly conducting 
inclusions Bj = B(cj,rj), j = 1,2. We assume that the inclusions are away from d£l, 
namely, there is a constant cq such that 

dist(Bj,dn) > c , j = 1, 2. (3.1) 
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We consider the following boundary value problem: 



' Au = inft\£iU B 2 , 
du 

"3 =9, 
ov an 

u = constant on dBj, j = 1,2, 



Ji ds = 0, j = 1, 2. 

dBi dv 



(3.2) 



Here g £ Lq(Q) (0 indicates that fg^g = 0) and we impose the condition that J 9 qU = 
for the uniqueness of the solution. 

In this section we derive an asymptotic formula similar to (|2.10p for the problem (|3.2D . 
Here we only consider the Neumann problem. But the same arguments work equally well 
for the Dirichlet problem. 

Let Aft : Lq(<90) — > H 1 (di}) be the Neumann to Dirichlet (NtD) map, i.e., 



where u is the solution to (|3.2p . Because of the assumption (|3.ip . we have 

ll A nb]||fl-i(9a) < C||5lU 2 (9n) 



(3.3) 



(3.4) 



for all g G Lq(8Q) for some constant C independent of e. See for example [13} Theorem 
2.2] 

For a bounded domain B with C 2 boundary let 5b and T>b denote the single and 
double layer potentials on B: 



S B [<p](x) = ^ / ln|x-y|^(y)ds(y), x e M 2 , 
(x-y,i/(y)) 



PbM(x) = - — 



¥>(y)cfo(y), xeK 2 \9B. 



We note 5b maps, as an operator defined on dB, C 0,a (dB) into C 1,a (dB) if a > 0. Thus if 
ip e C°' a (aB), then «S B [y>] belongs to C x ' a (B) and C 1 ' a (R 2 \ B). The single layer potential 
enjoys the following jump relation 



d(S B [<p\) d(S B [ip]) 



du 



Du 



(p on dB. 



(3.5) 



It is known that there are harmonic functions H and a pair of potentials {ipi,^) £ 
Cq'°(9-Bi) x Cg' a (9i?2) (0 indicates that the integral of <pj over dBj is zero) for some a > 
such that the solution u to (|2.10p is represented by 

n(x) = F(x) + 5 Bi [^ 1 ](x)+5b 2 ^ 2 ](x), x6fi\(B 1 UB 2 ). (3.6) 

In fact, H is given by 

tf(x) = -S n [9](x)+Pn[A n [g]](x), x e fi, (3.7) 
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and (if\, if2) is the unique solution to 

, d(S B2 [<p 2 ]) 8H 

d(S B M) , . dH 

dvw + X<p2 = ~mM ondB2 > 

where A = |. Here v^>> denotes the normal vector to j = 1,2. See [144 I15| (also 

mm)- 

Let Oi and O2 subdomains of O such that f^i C O2 and ^ C H. We further assume 
that B\ and B 2 are still away from d£l\, i.e., 

dist(J3i U £ 2 ,0ni) > ci, (3.9) 

for some ci > 0. By the Runge approximation, there is a sequence of harmonic functions 
H n in M 2 such that H n ^ H as n — > 00 in L°°(fi2)- For each n, let u n be the solution to 
(|2.4p with i7 replaced with Then u n can be represented as 

n„(x) = £ n (x)+cS Bl [^™ ) ](x)+cS i?2 ^ 2 n) ](x), x€l 3 \(BiUB 2 ), (3.10) 

where (</?i ,^4 ) * s unique solution to (|3.8p with H replaced with H n . Since — > 
as n —7- 00 in C°' a (dBj) for j = I and 2, we infer from the linearity of the integral 

equation that {(p^\<p^) -> (pi, y> 2 ) as n -> 00 in C°- a (aBi) x C°' a {dB 2 ). It means 
that u n — > u in C 1,a (Qi \ (B\ U £2))- We thus get the following theorem from Theorem 
I2TT1 

Theorem 3.1 Let u be the solution to (|3.2p and H be the function defined by (|3.7|) . Then, 
the solution u can be expressed as follows: 



u(x) = ah(x) + 6(x), xGfi\5iUB 2 (3.11) 

w/iere 

a= 4vrrir2 
n + r 2 

and 

l|V6|| Loo(OV(BlUB2)) < C (3.13) 
/or a constant C independent of e. 

It is worth looking more closely at the formula (|3.12p of the stress intensity factor. 
The function H is given by 

H{x) = ~[ ln|x-yb(y)d S (y)-i / (x ~ y '" ( 2 y)) A n [g](y) ds(y), 



and hence 



2-7rrir 2 
vr(ri + r 2 ) 



<P-y.n) / n , / x 
— f2-5(y) ds(y) 



en IP — y| 

- v nWr. - v „<V|\n 1 

(3.14) 



+/J^- (p " y ;V Ky) V ^>^> 

So if we can measure the Dirichlet data [g] on 90, we can determine the intensity of the 
stress using the boundary data. We emphasize that a is bounded regardless of e thanks 
to (EX 
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4 The insulated case 



We now deal with the case when circular inclusions are insulated, i.e., the conductivities 
are 0. Consider the solution to the free space problem: 



Au = 
du 
dv 







in R 2 \ B 1 U B 2 , 
on 8B 1 UdB 2 , 
as Ixl — > oo. 



(4.1) 



tu(x)-F(x) = 0(|xr 1 ) 
From the jump formula of the single layer potential, u can be represented as 
tt(x) = #(x) +S Bl [<pi](x) +SbM(x), xeR 2 , 



(4.2) 



for a pair of potentials (cpi, ip 2 ) G Lq(8Bi) x Lq(8B 2 ) satisfying (|3.8p with A = —\. 

Let H be an harmonic function in R 2 such that H is a harmonic conjugate of H. 
Then the solution u to (|4.1|) is a harmonic conjugate in M 2 \ B\ U -B 2 of u which is the 
solution to (|2.4p with in the place of H, see for example [3]. Note that by the Cauchy- 
Riemann equation, the tangential derivative of u is the same as the normal derivative of 
u on the disks, and hence u is constant on each disk Bj, j = 1.2. Theorem 12.11 yields . for 
x G R 2 \ BiUB 2 , as e -)• 0, 



Vu(x) 



2rir 2 
ri + r 2 



;n-VF)(p) 



x - pi 

|x - pi r 



X - p 2 

|x - P2I' 



+ 0(1). 



Let t is the unit vector perpendicular to n such that (n, t) is positively oriented and 



x 



-x 2 



for x G R 2 . Since Vn = (Vu) and n • Vif(p) = t • X7H(p), we have 



Vit(x) 



2nr 2 
n + r 2 



(t-v#)(p) 



(x-pi)- 1 (x-p 2 ) 



x - pi 



x - p 2 



+ 0(1). 



(4.3) 



Using (|4,3p we can obtain an expression of the solution n to (|4.ip . Let arg : ]R 2 \ 
{(0,0)} — )• [— 7r, 7r) be the argument function with a branch cut along the negative real 
axis, where x = (xi,x 2 ) is identified with x\ + ix 2 . Define 



/i_l(x) = — ( arg(x - pi) - arg(x - p 2 ) - arg(x - d) + arg(x - c 2 ] 



(4.4) 



where Cj is the center of Bj, j = 1,2. Note that h± is a harmonic function well defined 
in R 2 \ {B\ U B 2 ) since the jump discontinuity of the argument crossing the branch cut is 
canceled out owing to Pj,Cj G Bj. Moreover, we have 

V/i_l(x) = ^-V(arg(x - p x ) - arg(x - p 2 ) - arg(x - c x ) + arg(x - c 2 )) 

+ 0(1). 




x - Pi, 



P2)" 



x - Pi 



x-p 2 
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Similarly to the free space case, the solution u to the boundary value problem with 
the insulated inclusion becomes the solution of the perfectly conducting disk by taking its 
conjugate. To be more precise, if u is the solution 



Au = infi\BiU B 2 , 
du 

— = ov.dBi\JdB 2 , (4.5) 



— = g on dQ, 
< ov 

where Q is a simply connected bounded domain C 2 -boundary and g E Lq(8Q), Then we 
have (g^D and $£5$ with A = -\ and 

fT(x) = -5 n [^](x)+X)nMan](x), xe!1. 

Since O is a simply connected domain, the harmonic function H admits a conjugate 
function H in 17. Similarly to in free space, there is a harmonic conjugate u of u in 
I 2 \ 5i U -B2 satisfies (|3.2|) with a harmonic conjugate in the place of iif. 
Thus, we have the following theorem. 

Theorem 4.1 Lei u be either the solution to j4-l\ ) or the solution to f^.5| j in which case 
H is the function defined by 7\). Then, u can be expressed as follows: 



u(x) = aj_/ij_(x) + 6j_(x), x outside B\ U B 2 (4.6) 

where 

»"l + r 2 

and 

||V6±|U» ( n\ (fllUfla)) < C (4.8) 
for a constant C independent of e. 



5 Numerical computations 



In this section we compute numerically the solutions to (|2.4[) and (|4.1|) . Computation of 
the solution in the presence of closely located inclusions with conductivity k = or oo is 
known to be a hard problem. To understand this difficulty, let us consider a standard way 
of computing the solution using the boundary integral method. 
The solution to (j2.4[) can be represented as 



u(x) = #(x) + S Bl bi](x) +«Sb 2 [p 2 ](x), xel 2 \ (BiU5 2 ) 



(5.1) 



where (<pi,cp 2 ) is the solution to (|3.8[) with A = i. We can compute (921,1^2) numerically 
by discretizing (|3.8p with M number of equi-spaced points on each disks Bi, i = 1,2. Let 
3C, = 1, . . . , M, be the nodal points on (9£>j and set 



A 



Mm 


A 12 




Vi' 


A 21 


A/m. 


^2 



(5.2) 
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where 



8H 



8H 



§H 

"cM 2 ) 



a„(i)v x i / 

and A\2 and A21 are the evaluation of the kernel of 



8H /„M> 

"blM ^ 2 J 



(5.3) 



^tjSb 2 and q^S Bi at nodes on 
dB\ and $-B 2 , respectively. We then obtain 1,(^2) by solving 



A 



<P2 



Y. 



(5.4) 



As Figure [T] shows, the matrix A has small singular values, and the condition number 
of A becomes worse as e tends to 0. Moreover, derivative of the kernel of jj^SBi [y>i](x), 
which is of the form ^ > is as big as ^ if x and y are on the arcs of dB\ and dB2 

which are close to each other. Hence, if ipi takes large values on those arcs, the error in the 
discretization of the single layer potential becomes significant. From Theorem 12.11 <pt is 
as big as when n • VH 7^ at the middle point of the shortest line segment connecting 
dB\ and dB2- Therefore, we need finer grids as e gets smaller, see Figured! 

We will show that this difficulty can be overcome by using the characterization of 
singular terms given in (|2.10|) and (|4.6p . 
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Figure 1: The first graph shows the singular values of A in the decreasing order (n) when 
e is 0.0020. The second graph shows the condition numbers of A as e tends to (from 
right to left). We use 256 grid points on each Bi, i = 1, 2, and hence the dimension of A 
is 512 x 512. 



5.1 Computation for the perfectly conducting case 

In this subsection we present a new method of computing the solution to (|2.4p based on 
the characterization of the singular terms obtained in this paper. 
Let 

h(x) = h(x) - — (log|x- ci| - log|x - c 2 |) 

= -!-(log|x- pi| -log|x-p 2 | -log|x-ci| +log|x-c 2 |) (5.5) 
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for x G M 2 \(SiUi?2)- This modified function has the property that fg B . h = for i = 1,2, 
which is useful for the computation. 

In view of (|2.10p . we look for a solution in the following form 



u(x) = aft(x) + H (x) + S Bl (x) + 5 Ba ty> 2 ] (x) 



x G M 2 \ (B l UB 2 ), 



(5.6) 



instead of (|5.ip , where a is given by (|2.1ip . According to Theorem 12. 1\ the gradient of 
the function H + Sb 1 [4>i] + 5b 2 [^2] is bounded on \ (B\ U -B2) for some bounded set Q 
containing B\ U-B2, and hence H^i ||L°°(aBi) and ||V , 2||L° o (aB 2 ) are bounded regardless of e. 

To find the integral equation for density functions (^±,ip2), we argue as follows: Let 
h e be the extension of h defined by /i e (x) = /t(x) for x G M 2 \ {B\ U B2) and 



7 1 logn 1 , . 

Z7T Z7T 

h \dB 2 + —^. — log|x-Ci|, X G i>2- 

Z7T Z7T 



(5.7) 



Then h e = h on dBi for vi = 1,2, and h e is harmonic in B\ and -B2 as well as in M. 2 \B\ U B2. 
Define 

u e (x) = a/i e (x) + iT(x) +5 Bl M(x) +5 B2 [^ 2 ](x), x G M 2 . 

Then u e is continuous in R 2 and harmonic in B\ U B2 ■ Since u e = u is constant on dBi , 
i = 1, 2, u e is constant in Bi, i = 1, 2. By taking the interior normal derivative of it e , one 
can see that 001,^2) 6 LlidBx) x Ll(dB 2 ) is the unique solution to 



1 dOSssM) 

Wij) i, ( 
^( 2 ) + 2^ 2 



<9# 



dh e 
dh e 



on dBi, 
on 9i?2- 



(5.8) 



Moreover, we have from ([5 
dh e 



dh e 



1 (X-C 2 ,^)(X)) 

_( x ) = — — - — — — 1 xedBi, 



2vr 



c 2 



.(*) 



1 (x- Cl ,^ 2 )(x)) 



2tt 



x G 0B 2 



X — Ci 



We can discretize (|5.8p and solve (j5.4|) to obtain (ipi,^)- Here Yi and I2 are given by 



8H 



8v( 



a (xj-c 2 , ^(xj)) 
2^ |xj-c 2 | 2 



9H_ U Mx _ a (xf-c 2 , ^W(xf )) 
1) V x l J 2tt 1-^- 



|xf-c 2 | = 



(4) + |f 



a (xj-d, ^ 2 >(x 2 )} 



|x 2 - Cl | = 



o <xf- C1) ^( 2 )(xf )} 



|x«-cip 



(5.9) 

While (<pi,<p2) m the representation (|5.ip increases arbitrarily as e tends to 0, (^1,^2) 
stays bounded. The difference between the actual (1^1,1^2) and the numerically obtained 
one is much smaller than that for ((^1,992) as Figure [3] shows. 
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The first graph of Figure [2] shows the inner products of Y in (|5.4p with singular vectors 
of A corresponding to small singular values. The dotted graph is when (|5.3|) is used and 
solid one is when (|5,9p is used. The inner product using (|5.3p is larger than the one 
using (|5.9p . This is expected: since the difference between the single layer potential and 
its discretization is large in the narrow region between two disks, the singular vector 
(corresponding to small singular values) components of this difference is not small. The 



second graph in Figure [2] shows the inner products of A 



y"2 



Y with singular vectors 



of A corresponding to small singular values, where (^,^2) is the (numerical) solution to 
(|5.4p using Y in (|5.3p (dotted line) and in (|5.9p (solid line). This numerical solution is 
obtained using the method described in subsection 15.31 (with high precision). The graphs 
clearly show that the method of this paper works much better than the standard method. 
Here H(x) = x\. 
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Figure 2: The first graph shows the inner product of singular vectors (corresponding small 
singular values) of A with Y in ()5.4p . where Y is the discretization of the right hand side 
of (j3.8[) (the dashed line), and of ()5.8[) (the solid line). The second graph shows the inner 
product of singular vectors (corresponding small singular values) of A with the error in 
(|5.4p . The radii of disks are fixed as T\ = T2 = 1, and the number of equi-spaced grid 
points is 256 on each dBi. The distance e = 0.0020, and the background potential is given 
by H(x) = x\. n indicates the location of singular values when listed in decreasing order. 



5.2 Computation for the insulated case 

Let h± be the function defined by (|4.4p . Since arg(x — pi) — arg(x — P2) — arg(x — ci) is 
a harmonic conjugate of log |x — pi | — log |x — P2 1 — log |x — ci | , which is constant on dB\ , 
we have 

d 

— (arg(x - pi) - arg(x - p 2 ) - arg(x - ci)) = 0, on dB ± . 



dv 



Hence, 



dh± 



x 



1 <9(arg(x - c 2 )) 



2tt 



1 <(x-c 2 )V (1) (x)> 



dBx 2ir 



x 



c 2\ 



x e 8B X . 



(5.10) 
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Similarly, we have 



dh± 



(x) 



1 <9(arg(x — ci)) 



1 <(x- Cl )V( 2 )(x)> 



8B 2 



2tt 



X — Ci 



dv( 2 ) v ' 2vr du( 2 ) 
We look for a solution u to (|4,ip in the following form: 

u(x) = a ± h ± (x) + fl"(x) + 5 Bl [^i](x) + 5 B2 [^ 2 ](x), x G M 2 \ (Bi U 5 2 ). 
where a± is given by (|4.7p and (tpi,ip2) £ Lq(8Bi) x Lf i {dB2) is the solution to 



x G dBa. (5.11) 



1 



d(<W^i]) 



8H a ((x-c 2 ) ± ,i/ (1) 



x)) 



2tt 



x 



c 2 



on <9.E>i, 



a ((x- Cl )^,^ 2 )(x)) 



<9z/( 2 ) 2vr 



on 



c ii 



(5.12) 



(5.13) 



5.3 Numerical Illustration 

In this subsection, we illustrate results of numerical computations using the algorithms 
proposed in the previous subsections. Two discs are Bj = B(cj,rj), j = 1, 2, of radius rj 
and centered at ci = (— n — e/2,0) and c 2 = (r 2 + e/2,0). 

We compute the solution in two different ways and compare them to demonstrate the 
effectiveness of the method proposed in this paper. We first compute the solution using the 
standard representation of the solution, namely, we use ()5.ip and solve numerically (|3.8p . 
The discretization for the computation was described at the beginning of this section. We 
denote by u the solution computed by this method. We then compute the solution using 
the representation (|5.6p and solve (|5.8p . The solution is denoted by u h . For comparison 
we solve (|5.8p yet another method which provides the solution with higher precision (but 
with high cost). 

Let Ri, i = 1,2, be the reflection with respect to the disks Bi defined by (|2.6p . We also 
define the the reflection of a function / by (i?^/)(x) = /(i?j(x)) for x G M 2 . Using the 
same argument as in [3] (see also [H]), one can show that the solution (^i,Y> 2 ) to (|5.8p is 
given by 



m=0 

oo 

-2 V 



8 



lj'2 







(i? 2 i?l 



H H log |x 

2vr &l 



c 2 



i2 2 [tf-Alog|x-ci|] 



m=0 



( Rl R 2 r(H-— log|x- Cl 

Z7T 



ill [H + — log|x-c 2 

Z7T 



8Bi 



8B 2 



We denote by u R the solution obtained by this method. We compare these solutions 
for various values of e. The radii are fixed as r\ = 2 and r 2 = 1.5, and the number 
of equi-spaced grid points is 256 on each dBi. The background potential is given by 
H(x) = 2x\ + (x 2 — x 2 ). Figure [3] shows the relative L 2 -errors of the normal flux ^ and 

jjjj- , compared to the normal flux . The vertical axis in the figure represents the values 
of 



du h 



du R 



L 2 (9Bi) 



2 


du R 








L*{dB x ) 



+ 



du h 



L2(8B 2 ) 



2 


du R 






dv<?) 


L*(dB 2 ) 
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indicated by circles and 



du du R 


L 2 {dB L ) , 


du du R 


L 2 (dB 2 ) 


2 


du R 


I 

i 2 (9Bi) 


2 


du R 
cM 2 ) 


L 2 (dB 2 ) 



indicated by squares. As e decreases (from right to left), the relative error increases for u 
and u h as expected. However, the relative error of u is notably small compared to that 
of u: when e ~ 0.001, the relative error of u h is as small as 0.01, but that of u is as big as 
1. 
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□ — 
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Figure 3: The circles are the relative errors of compared to and the squares 

are those of Both depend on the distance e between two perfectly conducting disks. 
The solution based on the asymptotic expansion has much smaller error. We use 256 grid 
points on each Bi, i = 1, 2. 

In Figure^ fixing e = 0.0156, we compare the relative errors for different grid numbers. 
Both radii are 1, and H(~k) = x\. The difference between the normal flux ^ and 
takes its maximum value at the point nearest the middle point of the shortest line segment 
between dB\ and dB2, which is the same for ^fjy- As the grid numbers increase, both 
the relative L 2 -error and the maximal difference decrease. But, relative errors of u h are 
much smaller than those of u. 

In Figure [5] and [6l the uniformly spaced contour level curves are drawn for the free 
space conducting and insulating case, respectively. The distance e = 0.0156 and the 
number of grid points on each disk is 256. The radii are t\ = t<i = \ except the lower-right 
figure where T2 = 2. The entire harmonic function H(pc) = x\ in the upper-left and the 
lower-right figure, and H(x) = X2 in the upper-right one, and -ff(x) = x\ — X2 in the 
lower-left one. 
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ah fi R 

Figure 4: The circles are the relative errors of compared to 7577, and the squares are 
those of ^7 in L 2 and L°° norms for various grid numbers. The solution based on the 
asymptotic expansion has much smaller error. 
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Figure 5: Level curves of the free space conducting case. The entire harmonic function 
H(x) = x\ in the upper-left and the lower-right figure, and H (x) = %2 in the upper-right 
one, and H(x) = x\ — %2 in the lower-left one. 
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